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Abstract 

Rather  than  delivering  conventional  munitions  through  the  airspace  of  uncooperative 
nations,  a  constellation  of  space-stored  weapons  could  potentially  target  any  point  on  the 
Earth  and  arrive  within  the  time  it  takes  to  de-orbit  and  re-enter  through  the  atmosphere. 
The  research  involves  applying  the  dynamics  of  atmospheric  re-entry  to  a  Common  Aero 
Vehicle  (CAV)  and  defining  a  ‘footprint’  of  attainable  touchdown  points.  The  footprint 
is  moved  forward  to  create  a  swath  representing  all  the  possible  touchdown  points  in  a  90 
minute  window.  A  nominal  constellation  of  CAVs  is  established  using  a  ‘streets  of 
coverage’  technique,  and  both  analytic  studies  and  numeric  genetic  algorithm  techniques 
are  used  to  modify  the  nominal  constellation.  A  minimum  number  of  CAVs  is  identified 
which  ensures  payload  delivery  to  an  area  of  interest  within  90  minutes. 
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OPTIMAL  CONSTELLATION  DESIGN  FOR  ORBITAL  MUNITIONS 


DELIVERY  SYSTEM 


I.  Introduction 


Background 

The  Common  Aero  Vehicle  (CAV)  is  a  lifting  body  capable  of  atmospheric  re¬ 
entry  (1 :29).  This  weapon  platform  could  be  deployed  on  air-launched  suborbital 
missiles,  ICBMs,  or  launched  into  low  Earth  orbit  via  conventional  boosters.  The  CAV 
is  envisioned  to  be  self-guiding  toward  its  target,  using  inertial  and  possibly  GPS 
navigation  in  concert  with  aerodynamic  controls.  When  placed  in  orbit  around  the  Earth, 
it  could  be  used  to  deliver  a  munitions  payload  to  any  location  within  its  re-entry 
footprint.  Eurthermore,  a  constellation  of  such  vehicles  could  give  100%  delivery 
coverage  over  any  desired  portion  of  the  Earth’s  surface. 

Problem  Statement 

In  response  to  a  query  by  the  National  Security  Space  Architect  (NSSA),  we  will 
attempt  to  quantify,  both  analytically  and  numerically,  the  minimum  number  of  CAVs 
required  to  fully  cover  a  given  portion  of  the  Earth.  Terrestrial  delivery  is  required  to 
occur  no  later  than  90  minutes  from  the  time  a  decision  is  made  to  strike  a  target. 
Coverage  may  be  any  band  of  latitude,  extending  from  0°  to  the  latitude  of  interest. 
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Research  Objectives/Focus 

The  research  involves  several  disciplines  and  will  determine  optimal  solutions  for 
constellations  of  CAVs,  dependent  upon  several  design  parameters.  An  exploration  of 
atmospheric  re-entry  is  necessary  to  determine  the  touchdown  footprint  of  a  single  CAV. 
Analytic  constellation  design  will  be  used  extensively  to  define  several  types  of  baseline 
constellations.  Numeric  genetic  algorithm  (GA)  techniques  will  be  used  to  search  for 
non-analytic  solutions.  Finally,  we  will  compare  the  results  of  both  techniques  and 
identify  the  most  efficient  types  of  constellations  to  use  in  this  application. 

Methodology 

While  some  of  the  research  involves  analytical  evaluation  of  CAV  constellations, 
a  great  deal  of  the  work  depends  upon  the  results  of  numeric  simulation.  Footprint  width, 
a  fundamental  quantity  used  in  the  analysis,  is  solely  determined  Ifom  numeric 
integration  of  the  CAV’s  equations  of  motion.  Additionally,  generation  of  Earth 
coverage  statistics  as  well  as  the  entire  GA  routine  is  numeric  in  nature.  All  of  these 
numeric  techniques  are  carried  out  using  MATLAB©  (10),  with  the  GA  routine  using  an 
add-on  software  package  Ifom  Optimal  Synthesis©  (11). 

Assumptions/Implications 

Since  this  work  represents  a  first  look  at  this  combining  atmospheric  re-entry  and 
constellation  design,  there  are  several  basic  assumptions  which  were  made  in  order  to 
reduce  the  computational  complexity  of  the  problem  and  obtain  a  first-order  solution. 
First,  the  Earth  is  assumed  to  be  spherical  and  non-rotating.  Second,  the  atmosphere  is 
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assumed  to  be  exponential.  Third,  gravity  is  assumed  to  be  eonstant  throughout  the 
CAV’s  trajectory.  Fourth,  we  assume  the  CAV  does  not  have  any  delta-v  capability  other 
than  that  required  to  de-orbit.  Finally,  the  CAV  is  not  placed  under  any  heating,  dynamic 
loading,  or  g-force  constraints.  Application  of  these  assumptions  leads  to  a  more 
conservative  design  than  might  be  possible  using  more  complex  techniques. 

Preview 

Analytic  results  point  to  a  polar  inclined  streets  of  coverage  (SOC)  constellation 
as  being  the  most  efficient  way  to  obtain  100%  coverage  for  high  latitudes.  However, 

GA  techniques  reveal  a  modified,  inclined  SOC  constellation  that,  when  investigated 
further,  can  be  obtained  analytically  and  provides  an  improvement  over  the  polar  SOC 
constellation  when  certain  coverage  requirements  are  imposed. 
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II,  Literature  Review 


Overview 

There  has  been  a  great  deal  of  research  in  the  disciplines  of  constellation  design, 
genetic  algorithm  (GA)  search  techniques,  and  atmospheric  re-entry.  In  some  cases,  GA 
techniques  have  been  applied  to  satellite  constellations  (2:169-77),  but  the  two  disciplines 
of  atmospheric  re-entry  and  constellation  design  have  generally  been  treated  separately. 

Relevant  Research 

Much  of  the  research  in  constellation  design  has  focused  on  minimizing  the 
number  of  communications  or  remote  sensing  satellites  required  to  continuously  cover  at 
least  some  portion  of  the  Earth  (3,  4:179-84,  5:31-64,  6,  7:1419-30).  These  works  all 
begin  with  the  direct  relationship  between  swath  width  and  satellite  altitude.  Satellites 
are  assumed  to  have  circular  footprints,  and  analysis  consists  of  examining  the  number  of 
orbit  planes  and  the  number  of  satellites  per  plane  as  variables  leading  to  the 
determination  of  swath  width.  Once  swath  width  is  obtained  and  the  constellation  is 
minimized,  the  required  altitude  can  be  directly  calculated.  Conversely,  constraints  on 
orbit  altitude  may  dictate  a  swath  width,  which  can  then  be  used  to  find  a  minimum 
number  of  satellites  that  yield  the  desired  coverage. 

Constellation  design  using  the  streets  of  coverage  (SOC)  approach  has  also  been 
investigated.  In  this  method,  the  ascending  nodes  of  orbit  planes  are  evenly  spaced 
through  360°  for  arbitrarily  inclined  constellations,  and  through  180°  for  polar  inclined 
constellations  (7:1420,  8:188-200,  9:431-33).  These  works  arrange  satellites  such  that 
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their  circular  footprints  are  aligned  in  such  a  way  as  to  minimize  the  number  or  orbit 
planes  required.  This  is  generally  done  by  placing  the  ‘dip’  created  by  adjacent  circular 
footprints  next  to  the  ‘bulge’  created  by  a  footprint  in  the  next  orbit  plane.  However, 
because  the  last  orbit  plane  is  counter-rotating  with  respect  to  its  next  neighbor,  its  nodal 
spacing  must  be  smaller  than  the  average  spacing  (7:1420-23).  In  any  case,  none  of  these 
works  investigates  nodal  spacing  for  values  other  than  180°  or  360°. 

Some  research  has  also  focused  on  determining  the  intersection  of  both  co¬ 
rotating  and  counter-rotating  swaths  in  order  to  facilitate  coverage  of  specific  latitudes  or 
latitude  bands  (3:8,  5:62-3).  An  analytic  method  of  determining  the  latitude  of  swath 
crossings  is  developed  using  spherical  geometry.  We  refer  to  this  analysis  extensively  in 
the  analytic  portion  of  this  work. 

Genetic  algorithm  search  techniques  are  widely  researched  and  documented.  No 
new  techniques  are  presented  here;  rather,  standard  GA  search  techniques  (12:21 1-15)  are 
employed  with  the  aid  of  a  MATLAB©  (10)  add-on  software  package  from  Optimal 
Synthesis©  (11).  We  refer  the  reader  to  texts  by  Holland  (17)  and  Koza  (18)  for  more 
information  on  genetic  algorithms  in  general. 

Atmospheric  re-entry  is  also  a  well-researched  subject.  Many  theoretical  and 
practical  studies  have  been  conducted  on  hypersonic  re-entry  vehicle  dynamics  and 
control.  Generally,  these  studies  have  focused  on  recovering  manned  spacecraft  (13:239- 
68)  or  on  ballistic  re-entry  of  ICBM  warheads  (14:8-16).  Of  specific  interest,  however,  is 
a  controllable  re-entry  vehicle’s  footprint  of  possible  touchdown  points.  This  topic  has 
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been  addressed,  and  the  footprint  has  been  analytieally  determined  (15:207-10).  The 
results  of  this  particular  work  are  critical  to  the  problem  addressed  in  this  study. 

Applicability  of  Current  Research 

This  research  was  sponsored  by  the  National  Security  Space  Architect  (NSSA)  in 
response  to  a  query  regarding  potential  offensive  space  architectures.  The  solutions 
presented  in  this  paper  represent  a  first  look  at  this  problem  Ifom  the  standpoint  of  storing 
munitions  on-orbit. 

This  problem  differs  Ifom  previous  research  in  that  swath  width  is  no  longer  a 
function  of  altitude  or  the  number  of  satellites  per  plane,  but  rather  a  fixed  value 
determined  solely  by  the  re-entry  performance  of  the  CAV.  Furthermore,  the  system 
does  not  operate  instantaneously  as  with  remote  sensing  or  communications  platforms. 
CAVs  cannot  deliver  their  payloads  until  they  have  physically  passed  through  the 
atmosphere,  which  consumes  a  finite  amount  of  time. 

Therefore,  there  is  a  specific  requirement  on  the  time  until  delivery.  This  allows 
us  to  account  for  both  re-entry  time  and  spacing  of  the  CAVs  within  an  orbit  plane.  In 
this  problem,  delivery  must  be  within  90  minutes  from  the  time  of  de-orbit. 

We  will  also  investigate  SOC  constellations  in  which  nodal  spacing  takes  on  some 
value  between  180°  and  360°.  This  approach,  combined  with  the  non-continuous 
coverage,  creates  a  unique  problem  to  solve. 
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Ill,  Methodology 


Overview 

We  begin  by  simulating  the  equations  of  motion  for  the  CAV  during  atmospheric 
re-entry  to  obtain  a  maximum  lateral  distance  and  the  time  to  attain  this  distance.  With 
this  information  we  define  the  area  that  the  munitions  could  impact  within  90  minutes. 
The  remainder  of  the  problem  consists  of  arranging  a  constellation  of  CAVs  such  that 
their  touchdown  swaths  completely  cover  the  Earth. 

Much  of  the  work  was  numerical  in  nature,  and  several  functions  were  created  by 
the  author  to  aid  in  processing  data.  They  include  a  re-entry  profile  function,  which 
simulates  the  equations  of  motion  and  outputs  latitude  and  longitude  as  a  function  of 
time;  a  constellation  development  function,  which  creates  a  nominal  constellation  of 
CAVs  based  on  inputs  such  as  swath  width,  inclination,  and  desired  latitude  coverage; 
and  an  Earth  grid  function  which  is  used  to  calculate  Earth  coverage  statistics.  The  GA 
portion  of  the  analysis  relied  heavily  on  a  fitness  function,  which  incorporates  the  number 
of  CAVs  in  a  constellation  and  the  percentage  of  Earth  coverage  generated  by  the 
constellation  to  produce  a  fitness  value  relative  to  all  other  constellations  being 
considered.  Einally,  the  GA  algorithm  used  many  built-in  functions  included  in  a 
MATEAB©  add-on  package  Ifom  Optimal  Synthesis©.  See  Appendix  C  for  the 
MATEAB©  code  used  in  these  functions. 
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Atmospheric  Re-entry 

We  begin  by  defining  the  referenee  frame  in  whieh  the  CAV  operates.  Starting 
from  an  inertial  frame  X-Y-Z,  with  its  origin  at  Earth  eenter,  we  introduee  a  rotating 
frame  x-y-z,  also  with  its  origin  at  Earth  eenter.  This  frame  is  rotated  through  two  angles: 
Earth  east  longitude,  0,  and  Earth  latitude,  (p.  The  CAV’s  position  veetor  lies  along  the  x- 
axis.  A  third  frame  a-b-e,  also  rotating,  is  eentered  on  the  CAV.  The  a-b  plane  lies  in  the 
loeal  vertieal  plane,  with  the  b-axis  direeted  out  the  front  of  the  CAV.  The  e-axis  is  given 
by  e  =  a  X  b.  Erom  this  frame  we  define  the  flight  path  angle,  y,  measured  downward 
from  the  loeal  horizontal  to  the  veloeity  veetor;  the  heading  angle,  'P,  measured  from  the 
loeal  latitude  to  the  projeetion  of  the  veloeity  veetor  onto  the  loeal  horizontal;  and  the 
bank  angle,  a,  measured  from  the  a-b  plane  to  the  lift  veetor.  We  also  note  that  the  lift 
veetor  is  always  perpendieular  to  the  veloeity  veetor.  Eigure  1  shows  these  relationships. 
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Figure  1  CAV  Coordinate  Systems 


We  make  some  simplifying  assumptions  before  proceeding.  The  Earth  is 
assumed  to  be  spherical  and  non-rotating.  Additionally,  the  atmosphere  is  assumed  to  be 
exponential,  and  gravity  is  assumed  constant  throughout  the  trajectory.  Based  on  these 
assumptions  and  reference  frames,  the  equations  of  motion  for  atmospheric  re-entry  are 
as  follows  (16): 
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We  refer  the  reader  to  the  seetion  on  notation  for  explanation  of  these  variables. 

The  lift  veetor,  as  the  shaping  foree  of  the  re-entry  trajeetory,  is  eontrolled  by  the 
bank  angle,  a.  This  is  a  similar  approaeh  to  that  used  in  the  Spaee  Shuttle  program  (13). 
Starting  from  a  point  immediately  after  the  re-entry  burn,  a  footprint  of  possible  impaet 
points  is  eonstrueted.  The  maximum  downrange  eapability  is  obtained  by  maximizing 
lift  and  holding  bank  angle  eonstant  at  0°.  Lateral  range  is  obtained  by  eommanding 
bank  angle  to  some  value  other  than  0°  in  order  to  give  the  lift  veetor  a  horizontal 
eomponent,  whieh  then  turns  the  vehiele  through  its  deseent.  Optimal  eontrol  of  bank 
angle  in  maximizing  lateral  range  has  been  investigated  (15:208),  but  for  this  effort  we 
ehoose  a  eonstant  bank  angle  to  simplify  the  proeess.  Vinh  gives  45°  as  a  suboptimal 
eonstant  value,  but  also  notes  that  for  any  given  lift-to-drag  ratio,  a  value  greater  than  45° 
will  produee  the  greatest  lateral  range.  This  optimal  value  is  obtained  by  solving  the 
eubie  equation 
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where  E  =  Cd/Ci  and  a  =  eos^  cr  (16:353). 

Onee  bank  angle  is  obtained,  the  equations  of  motion  are  numerieally  integrated 
in  MATLAB©  using  a  variable  step  size  5^*'  order  Runge-Kutta  algorithm.  Of  particular 
interest  in  this  application  is  the  CAV’s  crossrange,  or  lateral  capability.  The  maximum 
lateral  range,  Xmax,  is  primarily  a  function  of  bank  angle  and  vehicle  ballistic  parameters. 

Although  control  in  this  simulation  is  open  loop,  we  choose  to  employ  a  simple 
method  of  control  to  help  maximize  lateral  range.  In  this  scheme  the  CAV  maintains  its 
optimum  bank  angle  (from  Equation  2)  until  the  heading  angle  is  turned  90°  away  from 
the  initial  heading,  at  which  time  the  bank  angle  is  set  to  0°.  This  prevents  the  CAV’s 
trajectory  from  becoming  a  spiral  and  allows  a  greater  lateral  range  than  if  the  bank  angle 
were  fixed  throughout  the  trajectory.  We  also  obtain  the  time  to  reach  Xmax,  denoted  as  fre- 
entry,  and  the  downrange  distance  oiXmax,  denoted  as  dre-entry- 
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Figure  2  Footprint  Size  vs.  Lift-to-Drag  Ratio 

There  are  several  quantities  from  Equation  1  which  might  change  the  outcome  of 
the  simulation.  Chief  among  these  are  the  ballistic  quantities  associated  with  the  CAV 
itself;  its  mass,  frontal  surface  area,  and  its  coefficients  of  lift  and  drag.  We  assume  a 
fixed  mass  and  frontal  area,  and  focus  instead  on  the  effects  of  the  ratio  of  lift  to  drag, 
denoted  as  L/D.  Figure  2  illustrates  the  effect  of  L/D  on  footprint  size  and  displacement 
from  a  de-orbit  burn  at  0°  latitude  and  0°  longitude. 

The  left  ends  of  the  footprints  are  open  due  to  the  fact  that  we  do  not  allow  bank 
angle  to  change  with  time  (e.g.  performing  roll  reversals).  The  complete  footprint  can  be 
obtained  using  more  sophisticated  methods  (15:207-10),  but  for  this  application  we  are 
only  interested  in  the  maximum  width  of  the  footprint.  Table  1  lists  some  possible  values 
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for  L/D  and  the  associated  values  for  lateral  range,  Xmax,  time  of  re-entry,  Ue-entry,  and 

downrange  distance  of  Xmax,  dre-entry- 


Table  1  L/D  and  Swath  Dimensions 


L/D 

^max 

tre-entry  (SCC) 

dre-entry 

0.3 

1.52° 

1921 

94.16° 

0.5 

3.75° 

2186 

103.77° 

0.7 

7.00° 

2518 

115.94° 

0.9 

11.25° 

2907 

129.91° 

Analytic  Constellation  Design 

We  begin  the  analysis  by  defining  the  total  area  which  can  be  covered  by  a  single 
CAV  within  the  90  minute  time  constraint.  Since  our  starting  point  can  be  anywhere  in 
the  CAV’s  orbit,  we  define  a  swath  of  coverage  based  on  the  current  position  of  the  CAV 
within  its  orbit,  wo;  the  orbital  mean  motion,  co;  and  the  quantities  Xmax,  t re-entry,  and  dre-entry- 

The  swath  length  is  defined  as  follows.  The  CAV  can  attain  any  point  within  the 
footprint,  but  we  are  interested  only  in  the  points  at  which  the  swath  is  at  its  widest. 
Therefore,  we  do  not  consider  points  before  or  after  dre-entry-  This  omission  gives  a 
conservative  estimate  of  the  ground  swath  (as  discussed  below),  but  greatly  simplifies  the 
analysis.  The  closest  point  along  the  ground  trace  is  given  by  uo  +  dre-entry-  If  we  allow 
the  CAV  to  travel  through  its  orbit  until  the  last  possible  moment,  defined  by  90  -  Ue-entry, 
we  obtain  the  furthest  point  along  the  ground  trace  that  the  CAV  can  attain.  The  swath 
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consists  of  the  area  between  these  two  endpoints.  Figure  3  illustrates  the  relationship 
between  time  and  distance  and  the  CAV’s  swath  length. 


Nominal  Orbit 
Trajectory 


Figure  3  CAV  Swath  Length 


Although  the  CAV  is  capable  of  attaining  any  point  in  the  footprint,  the  time  it 
takes  to  travel  to  that  point  is  variable.  To  simplify  the  analysis,  we  choose  a  constant 
value  for  Ue-entry  Figure  4  illustrates  the  time  difference  between  two  points  in  an 
example  footprint,  with  L/D  at  0.7.  The  difference  between  the  banked  trajectory,  which 
takes  2518  seconds,  and  the  straight  trajectory,  which  takes  2562  seconds,  represents  only 
a  1.7%  deviation.  The  time  difference  grows  with  increased  L/D;  the  difference  is 
approximately  9%  at  an  L/D  of  1 .2.  We  will  always  use  the  longer  time  to  represent 
tre-entry  in  Order  to  maintain  a  conservative  design. 
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Figure  4  Time  Difference  for  Straight  vs.  Banked  Trajectory 

The  ends  of  the  swath  are  somewhat  irregular  in  shape  due  to  the  possibility  of 
using  varying  amounts  of  bank  angle  during  the  descent  (see  Figure  2).  In  this 
application,  however,  we  assume  that  the  swath  is  of  constant  width  w,  where  w  =  . 

Additionally,  we  note  that  the  length  of  the  swath  is  given  by 

l  =  (3) 

Although  this  approach  does  not  maximize  the  full  potential  of  the  CAV,  it  allows  for  a 
simpler  analysis  of  constellation  coverage.  Figure  5  and  Figure  6  illustrate  this  concept. 
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Now,  based  on  the  length  of  a  swath  of  eoverage,  we  ean  direetly  ealeulate  the 


number  of  CAVs  required  per  orbit  plane  by 


s  =  ceiling 

V 


In 

~T 


y 


(4) 


where  /  is  given  in  Equation  3  and  ceiling  is  a  funetion  that  rounds  up  to  the  nearest 
integer.  Sinee  s  must  be  an  integer,  there  will  likely  be  some  level  of  overlap  between  the 
ends  of  the  individual  swaths  within  the  orbit  plane. 


Figure  5  Aetual  vs.  Simulated  Footprint  Comparison 
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Figure  6  Actual  vs.  Simulated  Swath  Length  Comparison 


Streets  of  Coverage  Constellation  Design 

Given  that  we  have  a  continuous  swath  of  coverage  of  width  w  for  each  orbit 
plane,  our  task  is  to  arrange  the  planes  such  that  we  cover  the  desired  portion  of  the  globe 
in  the  most  efficient  fashion.  To  begin,  we  adopt  a  streets  of  coverage  (SOC)  approach, 
in  which  we  ensure  equatorial  coverage  by  setting  adjacent  orbit  planes  close  enough  so 
that  they  leave  no  gaps  at  the  equator  (8:191-93,  9:431-33).  SOC  constellations  may  be 
arbitrarily  inclined,  in  which  case  the  ascending  nodes  are  equally  spaced  through  360°. 
They  may  also  be  polar  inclined,  in  which  case  the  ascending  nodes  are  equally  spaced 
through  180°  (9:431).  We  will  investigate  both  these  options  as  well  as  a  third,  modified 
type  of  SOC  constellation  in  which  the  ascending  nodes  are  equally  spaced  through  some 
value  between  180°  and  360°. 
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For  inclined  SOC  constellations,  we  are  free  to  ehoose  any  inclination  for  the 
orbit  planes  as  long  as  the  required  latitudes  remain  covered.  We  also  note  that  as  the 
orbit  planes  beeome  more  inelined,  the  swath  will  cover  a  larger  portion  of  the  equator, 
thereby  reducing  the  total  number  of  planes  required  to  eover  the  entire  equator.  Figure  7 
illustrates  this  concept  and  shows  the  swath  as  it  crosses  the  equator. 


Application  of  spherical  trigonometry  gives  the  erossing  width  by 


c  =  sm 


sm  w 


V  sm  z  J 


(5) 


The  number  of  orbit  planes  required  to  produce  an  inclined  SOC  constellation  is  given  by 
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p  =  ceiling 


(6) 


In 
c 

Of  course,  p  must  also  be  an  integer  so  we  round  up  and  aeeept  any  overlap  that 
oeeurs  between  adjaeent  planes.  At  this  point  we  have  defined  an  inelined  SOC 
eonstellation. 

Sinee  the  size  of  the  swath  is  direetly  related  to  the  L/D  of  the  CAV,  and  the 
number  of  planes  is  direetly  related  to  both  inelination  and  L/D,  it  follows  that  different 
eombinations  of  inelination  and  L/D  will  require  different  numbers  of  CAVs  for  inelined 
SOC  eonstellations.  Generally,  as  L/D  inereases,  the  number  of  planes  deereases;  and  as 
inelination  inereases,  the  number  of  planes  also  inereases.  Figure  8  shows  the  number  of 
orbit  planes  for  inelined  eonstellations  as  a  funetion  of  L/D  and  inelination.  Figure  9 
shows  an  inelined  SOC  eonstellation. 
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Figure  8  Number  of  Orbit  Planes  Required  for  Inclined  SOC 

Constellation 


Figure  9  Inclined  SOC  Constellation  (i  =  60°,  Xmax  =  10°,  p  =  18) 


For  polar  SOC  constellations,  Equation  6  may  be  simplified.  This  is  due  to  the 
fact  that  at  inclinations  of  90°,  the  aseending  and  descending  paths  of  the  CAVs 
completely  cover  the  globe.  Thus,  spaeing  the  orbit  planes  around  the  entire 
circumference  of  the  Earth  would  result  in  two  CAVs  traveling  opposite  directions  over 
the  same  ground  traee.  To  eliminate  this  redundancy,  we  distribute  the  orbit  planes 
around  only  half  the  Earth.  The  number  of  orbit  planes  is  given  by 


p  =  ceiling 


n 

w 


(7) 


In  this  case  we  see  that  as  E/D  inereases,  the  number  of  planes  decreases.  Eigure  10 
shows  the  number  of  planes  for  polar  eonstellations  as  a  function  of  E/D.  Eigure  1 1 
shows  a  polar  SOC  constellation. 
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Figure  10  L/D  vs  Number  of  Planes  for  Polar  SOC 


Figure  1 1  Polar  SOC  Constellation  (i  =  90°,  Xmax  =  10°,  p  =  9) 


Modified  Streets  of  Coverage  Design 


In  Figure  9,  we  note  that  although  we  do  indeed  have  eomplete  eoverage  in  the 
area  of  interest,  we  also  have  a  great  deal  of  redundant  eoverage.  In  the  interest  of 
economy,  we  might  consider  ways  to  eliminate  one  or  more  orbit  planes  from  the 
inclined  SOC  constellation  to  produce  a  modified  SOC  constellation. 

We  first  consider  removing  half  the  orbit  planes,  as  shown  in  Figure  12.  That  is, 
the  ascending  nodes  are  distributed  around  180°,  just  as  in  the  polar  SOC  constellation. 
There  are  now  large  areas  of  non-coverage;  in  fact,  the  only  latitude  fully  covered  is  the 
equator.  This  is  obviously  not  an  effective  solution,  so  we  next  consider  removing  a 
smaller  number  of  orbit  planes.  In  Figure  13,  only  three  orbit  planes  have  been  removed. 
Coverage  is  only  slightly  reduced  and  full  coverage  still  exists  nearly  to  the  inclination  of 
the  constellation.  This  result  is  promising,  and  we  now  consider  how  altering  the 
inclination  of  the  modified  constellation  affects  the  minimum  number  of  orbit  planes 
required  for  full  coverage  below  a  specified  latitude. 
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Figure  12  Modified  Inelined  SOC  Constellation  with  coverage  at 
0°  (i  =  60°,  2„«,=  10°,p  =  9) 


0  30  60  90  120  150  180  210  240  270  300  330  360 

Longitude 


Figure  13  Modified  Inclined  SOC  Constellation  With  Coverage  at 
~±52°  (i  =  60°,  10°,p=  15) 


We  first  investigate  how  far  apart  two  planes  ean  be  while  still  eovering  a  given 
latitude.  To  determine  this  value,  we  must  understand  how  the  spaeing  between  orbit 
planes  relates  to  the  intersection  of  their  swaths  (5:62-3).  Figure  14  illustrates  the 
geometry  involved. 


Side  View 


Figure  14  Swath  Intersection  Geometry 


Several  formulas  Ifom  spherical  trigonometry  may  be  applied: 


sm 


v2y 

sin«  = 


sinx  = 


tan  J 
tan/ 
cos/ 
cosd 
sin /I 
sin« 


(8) 
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We  also  know  that 


(j  +  x)  =  (j) 


req 


(9) 


where  cpreq  is  the  latitude  of  required  eoverage.  These  relationships  ean  be  combined  to 
eliminate  a  and  x  and  produce  an  expression  for  n: 


sin 

V 


n 

2 


y 


sin^zi^^^  cost -sin /I 
cos^^^^  sinz 


(10) 


The  complete  derivation  of  this  formula  is  given  in  Appendix  B. 

For  this  application,  we  define  cpreq  and  then  find  a  value  for  inclination  such  that 
the  number  of  orbit  planes  is  minimized.  An  inclined  SOC  constellation  in  which 
inclination  is  equal  to  the  latitude  of  interest  serves  as  a  baseline  Ifom  which  we  hope  to 
improve.  We  will  show  that  it  is  more  efficient  to  use  orbit  inclinations  that  are 
somewhat  higher  than  cpreq. 

To  create  modified  SOC  constellations,  we  must  find  the  smallest  value  of  n  such 
that  Equation  9  is  satisfied.  Stated  another  way,  we  are  seeking  a  nodal  spacing  such  that 
the  intersection  of  the  upper  boundaries  of  two  swaths  occurs  at  (preq,  as  shown  by  Point  O 
in  Figure  14.  We  choose  a  value  for  i  and  find  n  using  Equation  10.  The  range  over 
which  the  ascending  nodes  are  distributed  must  be  at  least  equal  to  n  but  less  than  In, 
meaning  the  additional  range  n  must  be  between  0  and  %.  This  value  n  is  then  inserted 
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into  Equation  6  to  determine  the  number  of  planes  required  for  full  eoverage  at  the 
desired  latitude  and  inelination: 


p  =  ceiling 


n  +  n 
c 


(11) 


The  numerator  in  Equation  1 1  refleets  a  nodal  eoverage  ofn  +  n  radians  rather  than  2n 
radians.  This  proeess  is  repeated  for  all  values  of  i  between  cpreq  and  90°  to  produee  a 
solution  curve  unique  to  this  particular  value  of  cpreq- 

The  results  of  this  process,  compared  to  inclined  SOC  and  polar  SOC 
constellations,  are  shown  in  Eigure  15  and  Eigure  16  for  two  different  swath  widths.  In 
these  examples,  the  curve  on  the  far  left  represents  the  number  of  orbit  planes  required  to 
create  inclined  SOC  constellations  over  the  full  range  of  possible  latitudes.  Each  marker 
along  this  curve  indicates  the  number  of  planes  required  at  a  specific  cpreq.  The  curves  to 
the  right  show  the  number  of  planes  and  inclinations  required  for  a  modified  SOC 
constellation  covering  that  latitude.  Each  of  these  is  marked  with  the  corresponding 
symbol  for  its  particular  (preq.  The  horizontal  line  represents  the  number  of  orbit  planes 
required  to  create  a  polar  SOC  constellation. 
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Figure  15  Orbit  Planes  Required  at  Various  Latitudes  for  Swath 

Width  =  7° 


Figure  16  Orbit  Planes  Required  at  Various  Latitudes  for  Swath 

Width  =  14° 


Numeric  Constellation  Design 

We  now  employ  Genetic  Algorithm  (GA)  techniques  in  an  attempt  to  further 
optimize  the  constellation.  Briefly,  GA  allows  for  an  examination  of  a  large  search  space 
using  techniques  borrowed  from  the  biological  processes  of  evolution  (12:207-10). 
Individual  variables  are  designated  as  genes,  and  the  genes  are  arranged  into 
chromosomes.  Each  chromosome  represents  one  particular  arrangement  of  the  problem 
variables  (in  this  case  representing  a  constellation  of  CAVs)  which  may  then  be 
manipulated  by  genetic  operations  such  as  mutation,  crossover,  and  inversion.  These 
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processes  are  pseudo-random  and  given  enough  time  and  appropriate  rules  for 
determining  the  fitness  of  eaeh  ehromosome,  will  converge  to  some  optimal  solution  (or 
one  of  a  set  of  pareto-optimal  solutions)  (2:2). 

Defining  the  Problem 

The  quantity  to  be  minimized  is  the  total  number  of  CAVs  required  for  100% 
eoverage  of  the  latitudes  of  interest.  A  constellation  is  deseribed  by  both  the  number  of 
CAVs  it  eontains  and  the  coverage  it  provides.  These  two  parameters  are  combined 
within  a  fitness  funetion  which  defines  the  total  effectiveness  of  the  constellation. 

In  all  GA  applieations  the  problem  must  be  represented  as  a  chromosome  with 
distinct  parts  that  can  be  manipulated  by  the  GA  processes.  Here,  we  choose  to  eneode 
the  constellation  as  a  fixed-length  binary  string.  Before  we  can  take  this  step,  several 
intermediate  encoding  steps  are  necessary  to  ensure  proper  operation  of  the  GA. 

Encoding  and  Decoding  the  Chromosome 

Each  CAV  in  the  constellation  is  fully  deseribed  by  its  orbital  elements:  semi¬ 
major  axis,  inclination,  eccentricity,  right  ascension  of  aseending  node,  argument  of 
perigee,  and  mean  anomaly.  In  this  application,  semi-major  axis  is  fixed.  Additionally, 
all  orbits  are  assumed  circular  which  means  argument  of  perigee  and  mean  anomaly 
become  undefined;  in  this  case  we  represent  the  CAV’s  position  within  the  orbit  using 
argument  of  latitude  (9:28-31).  Therefore,  eaeh  CAV  ean  be  fully  described  using  only 
inclination,  right  ascension,  and  argument  of  latitude,  and  any  constellation  of  CAVs  may 
be  defined  by  an  (A  x  3)  table  of  these  values. 
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One  approach  to  building  a  chromosome  from  this  table  is  to  simply  arrange  the 
rows  one  after  the  other  into  a  single  vector  of  length  3N.  However,  there  are  several 
shortfalls  with  this  method.  First,  when  performing  a  crossover  operation,  two 
chromosomes  must  be  cut  at  some  point  along  their  length.  The  chromosomes  then 
exchange  all  information  contained  after  the  cut  with  each  other.  To  retain  all  the 
information  for  each  CAV,  the  cut  must  occur  along  the  length  of  the  chromosome  in 
some  multiple  of  three.  Currently,  the  software  used  for  the  GA  process  is  unable  to 
enforce  this  condition,  and  as  a  result  any  offspring  from  the  crossover  operation  are 
likely  to  be  missing  some  information.  Rather  than  attempt  to  work  around  these  issues, 
we  adopt  a  table  lookup  approach. 

We  create  a  table  which  holds  all  possible  combinations  of  inclination,  right 
ascension,  and  argument  of  latitude  within  certain  fidelity.  Each  row  of  this  table 
represents  one  possible  set  of  values  describing  a  single  CAV.  Individual  genes  are  now 
reduced  to  a  single  integer  representing  the  appropriate  row  in  the  lookup  table.  Figure 
17  illustrates  the  difference  between  the  two  approaches. 
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Genetic  Structure  #1  (Real  Number  Values) 
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Each  gene  represents 
a  real  value  corresponding 
to  an  orbital  element  of  a  CAV. 
Three  genes  are  required  to  obtain 
complete  information  on  one  CAV. 


Genetic  Structure  #2  (Binary  Lookup  Values) 


01101001101001010001 

00000100001011011111 

00111111100010010000 

Each  gene  represents  a  row  in 
the  lookup  table.  One  gene 
contains  complete  information 
on  one  CAV. 


Figure  17  GA  Encoding  Scheme  Comparison 


We  construct  the  lookup  table  by  allowing  the  variables  of  interest  to  increment  in 
steps  of  .05  radians;  inclination  varies  through  a  range  of  n  radians  while  right  ascension 
and  argument  of  latitude  vary  through  27i  radians.  This  allows  for  a  thorough  search  of 
possible  configurations  while  avoiding  the  computational  overhead  of  an  extremely  large 
matrix.  The  size  of  this  table  is  given  by 


size  = 


'' 2n:\( 2;rY  n  ^ 


.05 


y.05  yf -05  y 


=  992,200 . 


(12) 


To  accommodate  this  many  values,  a  20-digit  binary  string,  which  may  take  on 
any  value  between  0  and  1,048,576  is  ideal.  Analytic  results  show  that  we  can  expect  to 
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never  have  more  than  100  total  CAVs  (using  current  L/D  values),  so  we  select  a  fixed- 
length  chromosome  of  20  x  100  =  2000  bits.  This  binary  approach  allows  easy 
manipulation  of  the  information  in  the  chromosome  by  the  genetic  processes. 

This  method  also  allows  the  algorithm  to  easily  remove  CAVs  Ifom  a 
constellation.  There  are  56,376  unused  rows  in  the  lookup  table,  and  their  contents  are 
assigned  to  the  null  set.  When  the  algorithm  references  one  of  these  ‘empty’  rows  during 
chromosome  manipulation,  the  interpretation  is  that  the  CAV  does  not  exist.  Likewise, 
removing  a  CAV  is  as  simple  as  setting  its  gene  to  a  value  that  references  an  ‘empty’  row 
in  the  lookup  table.  Since  the  goal  is  to  minimize  the  number  of  CAVs,  having  many 
copies  of  the  null  CAV  was  thought  to  be  desirable  in  increasing  the  likelihood  of 
removing  additional  CAVs  during  crossover  operations. 

Encoding  the  chromosome  now  consists  of  replacing  each  row  in  the  constellation 
matrix  with  its  corresponding  integer  row  number  Ifom  the  large  lookup  matrix,  then 
converting  each  integer  into  its  20-bit  binary  equivalent.  The  final  step  is  to  rearrange  the 
now  (100  X  20)  matrix  into  a  (1  x  2000)  row  vector  for  processing  by  the  GA  code.  To 
decode  processed  chromosomes,  we  simply  reverse  the  steps. 

Genetic  Processes 

Once  the  encoded  chromosomes  are  created,  they  are  processed  to  generate 
variability  in  the  design.  There  are  three  processes  which  occur  within  each  generation; 
selection,  modification,  and  decimation.  Selection  is  performed  by  choosing  one  or  two 
chromosomes,  either  at  random  or  in  proportion  to  their  fitness.  Modification  consists  of 
either  mutation,  where  each  bit  is  subject  to  inversion  (1  ^  0  or  0  ^  1)  based  on  a  set 
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probability;  or  crossover,  where  two  chromosomes  exchange  a  randomly  determined 
amount  of  genetie  material.  In  both  oases,  chromosome  length  is  fixed.  Decimation 
occurs  at  the  end  of  each  generation  and  is  used  to  eliminate  the  least  fit  members  Ifom 
the  population,  and  to  keep  the  population  size  at  a  manageable  level. 

Coverage  Evaluation 

A  straightforward  method  to  evaluate  oonstellation  coverage  is  to  distribute 
evenly  spaced  points  around  the  equator  and  around  each  line  of  latitude,  then  oheck  if 
eaoh  point  is  covered  at  any  point  in  the  simulation.  However,  if  the  same  number  of 
points  are  distributed  along  the  higher  latitudes  as  along  the  equator,  there  will  be  many 
more  points  in  the  polar  regions.  This  will  tend  to  skew  any  figures  of  merit  toward  polar 
eoverage,  which  may  not  be  desirable.  We  eliminate  this  problem  by  reducing  the 
number  of  points  along  any  line  of  latitude  by 

points  =  (points  Xcos  (f)  (13) 

We  also  inelude  a  random  starting  point  on  eaoh  line  of  latitude  to  prevent 
artificial  weighting  of  the  prime  meridian.  As  the  number  of  grid  points  increases,  so 
does  the  time  required  to  oompute  eoverage.  Tests  were  performed  with  grid  spaoing  as 
low  as  1°  and  did  not  show  any  appreoiable  inorease  in  acouraoy  over  larger  values  when 
determining  coverage.  A  grid  spaoing  of  5°  was  chosen  as  a  good  balance  between  grid 
fidelity  (1656  grid  points)  and  oomputational  effioiency. 
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Fitness  Evaluation 


To  evaluate  the  fitness  of  a  eonstellation,  eaeh  CAV  in  the  eonstellation  is 
propagated  through  its  swath.  Sinee  the  swath  is  defined  only  by  Imax,  grid  eoverage  at 
eaeh  time  step  ean  be  eheeked  with  a  simple  dot  produet  ealeulation  (see  Figure  6).  If  a 
grid  point  is  eovered  it  is  flagged,  and  after  all  CAVs  are  propagated  the  total  eoverage  is 
ealeulated.  The  eonstellation’ s  fitness  is  given  by 


/  = 


number  of  CAVs 
coverage  ^ 


(14) 


where  q  represents  a  variable  exponent  designed  to  penalize  ineomplete  eonstellations. 
This  fitness  value  is  passed  baek  to  the  GA  eode,  whieh  then  ranks  the  eonstellation,  and 
either  retains  the  eonstellation  for  future  generations  or  diseards  it  through  the  deeimation 
proeess. 

In  Equation  14,  the  integer  q  in  the  denominator  eontrols  the  rate  at  whieh  the  GA 
algorithm  eonverges  on  possible  solutions.  Constellations  with  less  than  100%  eoverage 
are  always  penalized  aeeording  to  the  amount  of  eoverage  they  have.  Early  in  the  seareh, 
we  keep  the  penalty  low,  and  therefore  q  is  small.  This  allows  a  larger  variety  of  genetie 
material  to  remain  in  the  population.  However,  as  the  seareh  proeeeds,  we  must  begin 
eliminating  those  eonstellations  with  less  than  100%  eoverage.  Thus,  we  inerease  the 
eoverage  penalty,  and  so  q  inereases  throughout  the  life  of  the  GA  proeess. 
Experimentation  with  the  algorithm  leads  us  to  set  ^  =  3  at  the  beginning  of  the  proeess 
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and  allow  it  to  increase  along  a  parabolie  eurve  until  the  end  of  the  proeess,  at  whieh  time 
q  =  20.  The  exponent  is  eomputed  by 


17 


{total  generations  - \f 


{current generation  -if  +3 


(15) 


where  the  value  of  total  generations  is  input  by  the  user. 
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IV.  Analysis  and  Results 


Chapter  Overview 

This  chapter  details  the  results  of  both  analytic  and  numeric  analysis  of  the  Earth 
eoverage  problem.  Minimal  constellations  are  diseussed  and  verified  in  both  methods, 
and  general  design  conclusions  are  made. 

Analysis 

To  assist  in  evaluating  the  performance  of  the  analytic  versus  the  numeric 
methods,  we  ehoose  the  following  CAV  properties:  L/D  =  0.7;  mass  =  1000  kg;  frontal 
surface  area  =10  m^.  All  CAV  orbits  are  circular  with  a  semi-major  axis  of  500  km. 
Bank  angle  was  ealculated  at  40.1°  using  Equation  2  .  The  following  values  were 
generated  by  the  simulation:  dre-entry  =  1 15.94°;  Ue-entry  =  2562  sec.;  Xmax  =  7°.  We  choose 
a  latitude  coverage  band  of  ±  65°  for  the  first  simulation,  and  ±25°  for  the  seeond. 

Analytic  Results 

We  now  develop  a  nominal  constellation  of  CAVs  based  on  these  values.  Swath 
length  equals  182.77°,  calculated  using  Equation  3.  Swath  width  is  given  by 
w  =  2(Xmax)  =  14°.  The  number  of  CAVs  per  plane,  5  =  2,  is  found  from  Equation  4. 
Number  of  orbit  planes  and  constellation  inclination  varies  widely  with  constellation 
type. 

Eor  the  polar  constellation,  p  =  14  from  Equation  7,  and  i  >  83°. 
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For  the  inclined  SOC  constellation,  p=  \2  and  i  =  25°  for  the  ±  25°  case;  and 
p  =  25  and  i  =  65°  for  the  ±65°  case.  The  number  of  planes  was  determined  from 
Equation  6,  and  the  inclination  was  set  equal  to  cpreq- 

For  the  modified  SOC  constellation,  p  =  %  and  i  =  19.5°  for  the  ±  25°  case;  and 
p=  \A  and  i  =  82.5°  for  the  ±65°  case.  The  number  of  planes  was  determined  using 
Equation  11,  and  the  inclination  was  determined  from  Equation  10. 

We  see  from  Figure  15  and  Figure  16  that  as  inclination  increases,  more  orbit 
planes  are  required  to  fully  cover  the  latitude  band  of  interest  when  using  an  inclined 
SOC  constellation.  This  is  because  the  width  of  the  swath  as  it  crosses  the  equator 
decreases  as  inclination  increases.  However,  for  the  modified  SOC  constellations,  a  more 
complicated  curve  results,  and  this  trend  is  actually  reversed  for  mid  to  high  values  of 
(preq.  Although  the  swath  width  at  the  equator  is  decreasing,  the  number  of  orbit  planes 
we  can  remove  increases,  driving  the  total  number  of  orbit  planes  down.  We  also  note 
the  impact  of  (preq  on  the  length  and  slope  of  each  curve.  As  expected,  high  latitudes  can 
only  be  serviced  by  high  inclination  orbits,  and  the  effect  of  increasing  inclination  is 
more  pronounced. 

Obviously,  removing  orbit  planes  from  an  inclined  SOC  constellation  improves 
the  efficiency  of  the  constellation.  For  lower-latitude  coverage,  these  modified  SOC 
constellations  are  the  most  efficient  way  to  cover  the  area  of  interest.  As  latitude 
increases  past  a  certain  value,  however,  polar  SOC  constellations  become  more  efficient. 
The  exact  latitude  where  this  transition  occurs  is  a  function  of  swath  width  and 
inclination  and  is  not  analytically  tractable.  With  sub-polar  latitude  coverage 
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requirements,  polar  SOC  eonstellations  generate  unneeded  eoverage  in  the  polar  regions, 
but  the  faet  that  they  require  only  half  the  orbit  planes  as  inelined  full  SOC  eonstellations 
make  up  for  this  relative  ineffieieney.  Therefore,  a  polar  SOC  eonstellation  is  the 
generally  the  best  method  for  providing  mid-  to  high-latitude  eoverage. 

If,  for  some  reason,  polar  orbits  are  not  aehievable,  the  modified  SOC 
eonstellation  still  provides  a  signifieant  advantage  over  its  inelined  SOC  eounterpart.  As 
shown  in  Figure  15  and  Figure  16,  eonstellation  effieieney  goes  down  as  inelination 
deereases,  but  still  remains  better  than  the  full  inelined  SOC  baseline. 

An  illustrative  example  is  to  eonsider,  Ifom  Figure  15,  the  eurve  representing  a 
(preq  value  of  35°.  In  this  ease,  the  full  SOC  solution  requires  approximately  29  orbit 
planes,  while  the  modified  SOC  solution  requires  approximately  25  orbit  planes,  at  an 
inelination  of  approximately  34°.  The  polar  SOC  solution  in  this  ease  requires 
approximately  26  orbit  planes. 

Numeric  Results 

The  GA  proeess  was  implemented  with  an  initial  population  eontaining  several 
inelined  SOC  eonstellations  for  inelinations  at  and  above  the  latitude  being  studied.  In 
order  to  provide  the  algorithm  enough  information  to  begin  a  valid  seareh  of  the  solution 
spaee,  approximately  100  randomly  generated  eonstellations  were  also  ineluded  in  the 
population.  Additionally,  an  ‘all  zeros’  ehromosome  was  added.  This  ehromosome 
represents  an  empty  eonstellation  (by  refereneing  the  null  rows  of  the  lookup  table),  and 
allows  the  algorithm  to  remove  CAVs  Ifom  a  eonstellation  via  the  erossover  operation. 
The  algorithm  was  run  multiple  times  using  a  variety  of  values  for  total  generations. 


39 


Tests  were  performed  using  up  to  10,000  generations  with  no  improvement  over  lower 
values.  A  value  of  1500  generations  was  chosen  to  give  a  good  balance  between 
computing  time  and  depth  of  search. 

In  each  GA  run,  the  two  best  constellations  and  the  two  worst  constellations  were 
plotted  for  analysis.  As  expected,  the  worst  constellations  were  random  assortments  of 
CAVs,  with  coverage  often  dropping  below  50%.  These  were  included  in  the  output  to 
verify  the  algorithm  was  keeping  a  large  variety  of  possible  solutions  in  the  population. 
The  best  constellations  were  similar  to  the  analytic  solutions. 

The  GA  algorithm  produced  interesting  results.  Numerically,  the  best 
constellations  found  were  similar  to  the  ones  obtained  analytically.  Figure  18  shows  a 
solution  at  freq  =  25°,  while  Figure  19  shows  a  solution  at  freq  =  65°.  The  freq  =  25°  case 
was  identical  to  the  analytic  solution,  but  in  the  (preq  =  65°  case,  the  GA  algorithm 
eliminated  one  additional  orbit  plane  from  the  analytical  solution,  which  lowered 
coverage  values  to  >  97%.  This  was  typical  of  the  GA  results  in  general;  in  most  cases, 
the  best  constellations  had  coverage  slightly  less  than  100%. 
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Analytic  vs.  Numeric  Comparison 

At  this  point  we  wish  to  summarize  the  results  of  both  analyses.  Although  the 
numerie  analysis  shows  that  further  reduetions  from  the  analytie  results  are  possible,  it 
does  not  represent  a  signifieant  savings.  The  general  result  is  that  polar  SOC 
eonstellations  are  the  most  effieient  way  to  eover  mid  to  high  latitudes,  while  modified 
SOC  eonstellations  beeome  optimal  at  lower  latitudes.  The  exaet  latitude  at  whieh  this 
oeeurs  is  mainly  a  funetion  of  swath  width  and  inelination  and  is  not  analytieally 
traetable.  Table  2  and  Table  3  eompare  the  results  of  the  analytie  and  numerie  analysis. 


Table  2  Constellation  Summary  for  Latitude  Requirement  of  25° 


Constellation 

Type 

Latitude 

Coverage 

Req’t 

Inelination 

Number 
of  Planes 

Number 

ofCAVs 

Coverage 

(%) 

Polar  SOC 

25° 

>83° 

14 

28 

100 

Full  SOC 

25° 

25° 

12 

24 

100 

Modified  SOC 

25° 

19.5° 

8 

16 

100 

GA  Low 

25° 

19.5° 

8 

16 

100 

Table  3  Constellation  Summary  for  Latitude  Requirement  of  65° 


Constellation 

Type 

Latitude 

Coverage 

Req’t 

Inelination 

Number 
of  Planes 

Number 

ofCAVs 

Coverage 

(%) 

Polar  SOC 

65° 

>83° 

14 

28 

100 

Full  SOC 

65° 

65° 

25 

50 

100 

Modified  SOC 

65° 

82.5° 

14 

28 

100 

GA  High 

65° 

83.1° 

13 

26 

97.1 

At  this  point  we  must  aeknowledge  that  the  numerie  results  are  somewhat 
disappointing.  This  may  be  due  to  the  laek  of  any  truly  new  or  interesting  solutions  or  it 
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may  be  due  to  shorteomings  in  the  GA  routine  itself;  without  further  investigation  we 
eannot  eonfirm  either  suspieion.  However,  the  results  ean  be  partially  explained  by 
addressing  some  eoneerns  regarding  the  operation  of  the  GA  routine.  There  were  several 
roadbloeks  to  overeome  in  this  formulation,  not  all  of  whieh  were  satisfaetorily  resolved. 

First,  the  fidelity  of  the  lookup  table  was  of  some  eoneern.  Although  eaeh  orbital 
element  was  ineremented  in  steps  of  .05  radians  (equivalent  to  2.86°),  a  finer  table  might 
lead  to  more  robust  results.  The  possibility  exists  that  an  optimal  solution  was  missed 
due  to  the  proper  value  not  existing  in  the  lookup  table.  However,  the  exeessive 
eomputation  time  assoeiated  with  a  larger  table  prohibited  its  use  in  this  study. 

Seeond,  ealeulation  of  the  eonstellation’s  fitness  value  was  of  some  eoneern.  The 
problem  of  properly  weighting  Earth  eoverage  was  addressed  early  in  the  researeh  by 
allowing  the  exponent  q  in  Equation  14  to  vary  throughout  the  simulation.  However,  the 
evolution  of  the  population  is  very  sensitive  to  the  value  of  ^  and  the  rate  at  whieh  q  is 
allowed  to  grow  during  the  simulation.  Without  further  experimentation  we  eannot 
definitively  state  that  the  fitness  funetion  ideally  ealeulates  the  true  fitness  of  the 
eonstellation. 

Einally,  the  issue  of  population  diversity  and  its  effeet  on  eonvergenee  of  the 
algorithm  must  be  addressed.  Originally,  the  initial  population  was  seeded  with  a  single 
inelined  SOC  solution,  an  ‘all  zeros’  ehromosome,  and  approximately  100  randomly 
generated  ehromosomes.  This  setup  produeed  eonstellations  remarkably  similar  to  the 
modified  SOC  eonstellations  diseussed  above.  It  was  reasoned  that  the  presenee  of  the 
‘zero  ehromosome’  allowed  the  algorithm  to  remove  individual  CAVs  in  the  inelined 
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SOC  constellation.  However,  when  a  different  mix  of  initial  chromosomes  was  used,  the 


routine  converged  to  a  eonstellation  containing  only  one  CAV.  This  dependence  on  a 
suitable  initial  population  is  a  documented  shortfall  of  GA  searches  in  general  (12:107). 
The  current  assortment  of  constellations  in  the  initial  population  was  generated  through 
extensive  experimentation;  there  may  be  a  better  mix  of  constellations  that  would  yield 
better  results. 
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V,  Conclusions  and  Recommendations 


Conclusions  of  Research 

For  mid  to  high  latitude  coverage  requirements,  polar  SOC  constellations  are  the 
most  efficient  method  of  providing  full  coverage  within  90  minutes  of  a  decision  to  de¬ 
orbit  the  CAV.  For  low  latitude  coverage  requirements,  or  in  circumstances  where  polar 
orbits  are  not  achievable,  modified  SOC  constellations  are  the  most  efficient  method. 

The  exact  latitude  at  which  this  transition  takes  place  cannot  be  obtained  analytically,  but 
is  mainly  a  function  of  swath  width  and  inclination.  Furthermore,  in  cases  where  less 
than  100%  coverage  is  acceptable,  additional  orbit  planes  can  be  removed  Ifom  the 
modified  SOC  constellation.  The  latitude  at  which  these  modified  SOC  constellations 
become  more  advantageous  than  a  polar  SOC  constellation  varies  with  the  swath  width  of 
the  CAVs  in  the  constellation. 

Significance  of  Research 

The  elimination  of  one  or  more  orbit  planes  Ifom  inclined  SOC  constellations 
results  in  launch  cost  savings  (fewer  launches  required)  as  well  as  overall  system  cost 
savings  (fewer  total  CAVs  required).  The  design  paradigm  addressed  here  is  valid  for 
constellations  which  do  not  require  continuous  coverage,  although  the  method  could  be 
extended  to  more  standard  applications.  If  a  ring  of  satellites  in  a  single  orbit  plane  can 
produce  a  continuous  swath  of  coverage  on  the  ground,  the  method  presented  here  may 
be  applied  to  design  of  the  constellation. 
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Also  worth  noting  is  the  realization  that  reducing  coverage  requirements  can 
eliminate  additional  orbit  planes  in  cases  where  the  cost  per  orbit  plane  far  exceeds  the 
value  of  complete  coverage.  A  trade  study  using  this  paradigm  could  be  conducted  in  the 
design  phase  of  almost  any  constellation  of  this  type. 

Recommendations  for  Future  Research 

In  the  future,  a  more  robust  study  of  CAV  re-entry  should  be  conducted  to  more 
accurately  and  completely  define  the  footprint  and  swath  size  of  the  vehicle. 

Optimization  of  the  re-entry  trajectory  could  yield  greater  swath  size  and  a  subsequent 
reduction  in  constellation  size.  Creating  an  algorithm  to  model  the  irregular  shape  of  the 
footprint,  and  thus  the  entire  swath,  would  further  maximize  the  potential  of  a  single 
CAV  and  lead  to  additional  reductions  in  constellation  size. 

Adding  a  delta-v  capability  while  on  orbit  would  also  change  the  CAV’s  footprint 
and  swath  size,  with  a  reduction  in  the  number  of  CAVs  required  being  a  likely  outcome. 
The  simplified  analysis  presented  here  is  not  sufficient  to  model  the  ability  of  the  CAV  to 
change  its  orbit  plane  before  re-entry.  Significant  further  study  is  necessary  to  include 
this  capability. 

The  fidelity  of  the  study  would  be  improved  by  adding  Earth  rotation,  J2,  and 
other  perturbations  into  the  model.  These  effects  will  not  change  the  overall  shape  of  the 
constellations,  but  will  provide  further  validation  of  the  concept  as  well  as  a  logical  link 
to  the  next  part  of  the  study. 

Follow-on  research  should  attempt  to  model  a  complete  system  of  CAVs  along 
with  their  timing  and  target  opportunities.  Some  specifics  would  entail  creating  an 
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algorithm  to  define  which  CAV  to  select  for  deployment  against  a  specific  target,  given 
the  current  time  and  time-on-target  information.  The  study  would  need  to  include  a 
robust  model  for  orbital  motion,  a  complete  description  of  re-entry  times  to  every  part  of 
the  footprint,  and  a  method  for  choosing  the  CAV  most  likely  to  arrive  at  the  target 
within  the  90  minute  time  constraint. 

Finally,  a  more  robust  and  reliable  GA  routine  should  be  implemented.  Many  of 
the  shortfalls  of  the  GA  routine  were  addressed  in  the  analysis  section  of  this  paper  and 
will  need  to  be  addressed  before  the  GA  search  can  be  considered  complete.  Although 
further  reductions  in  constellation  size  may  not  be  realistic,  this  endeavor  should  be 
undertaken  to  eliminate  any  doubt  on  the  matter. 
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Appendix  A 


Notation 


EOM  VARIABLES 


h  =  altitude 
6  =  Earth  longitude 
(j)  =  Earth  latitude 

V  =  Earth  —  relative  velocity 

Y  =  flight  path  angle,  measured  downward  from  local  horizontal 

=  heading  angle,  measured  from  local  latitude  to  the  projection  of  v 
onto  the  local  horizontal 

(7  =  bank  angle,  measured  from  local  vertical  to  the  lift  vector 


D 

m 


L 

m 


^Pm 


—h 

ll 


^Pm 


Pm 

^E 

H 

PO 

S 


m 

=  Earth  radius 
=  scale  height 

=  atmospheric  density  at  sea  level 

=  surface  area  of  reentry  vehicle  normal  to  velocity  vector 


OTHER  VARIABLES 
a  =  semi  —  major  axis 
0)  =  orbital  mean  motion 
i  =  inclination 

Q  =  right  ascension  of  ascending  node 
u  =  argument  of  latitude 
Amax  —  maximum  crossrange  capability 
w  =  swath  width 
I  =  swath  length 

c  =  swath  width  at  equator  crossing  for  given  inclination 
s  =  number  of  satellites  per  orbit  plane 
p  =  number  of  orbit  planes 
d  =  latitude  of  ground  trace  crossing 
X  =  latitude  of  upper  swath  boundaries  crossing 

n  =  ascending  node  spacing  between  orbit  plane  at  Q  =  180°  and  last  orbit  plane  used 
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Appendix  B 


Nodal  Spacing  Derivation 

We  are  seeking  an  analytie  expression  for  the  value  n  in  Equation  10  and 
Equation  1 1  which  is  used  to  determine  the  range  of  nodal  crossings  in  a  modified  SOC 
constellation.  This  derivation  is  based  on  the  work  of  Rider  (5:62).  In  Rider’s  work,  n  is 
known  and  k  is  the  quantity  being  sought.  Our  approach  differs  from  the  reference  that  k 
is  known  and  n  is  unknown;  and  that  n  is  defined  differently. 


Side  View 


We  reproduce  Eigure  14  here  and  begin  with  the  following  relations  from 
spherical  trigonometry: 
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sin 


^  n\  tan  d 
yl)  tan  i 


(Al) 


cos  i 

sin  a  = - 

cosd 


(A2) 


sin /I 

sin  X  = - 

sin  a 


(A3) 


We  also  make  use  of  the  trigonometrie  identity 


sin(a  +  h)  =  sin  a  eos  b  -  eos  a  sin  b 


(A4) 


and  note  from  Equation  9  that 


X  =  <f>,e,  -d 


(A5) 


Substitution  of  Equation  A5  into  Equation  A3  yields 


=  (A6) 

sm  a 


Inserting  Equation  A2  and  applying  the  identity  in  Equation  A4,  we  obtain 


sm  X  eos  d  .  ,  ,  ,  •  , 

- =  sm  (p  eos  d  -  eos  (p  sm  d 

eos  i 


(A7) 
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This  expression  ean  be  rearranged  to  produee 


tan  d  =  tan 


sin /I 

eos/eos^/i,,^ 


(A8) 


Substituting  Equation  A8  baek  into  Equation  A1  and  rearranging  terms  gives 


sin 

V 


n 

2 


y 


sin  eos  i  -  sin  X 
sinz 


(A9) 


This  relationship  gives  n  entirely  in  terms  of  the  known  values  cpreq,  k,  and  i. 
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Appendix  C 


MATLAB©  Code 


Re-entry  Simulation 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  REENTRY  SIMULATION.  CALLS  EOMS  TILL  LOR  INTEGRATION 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  [lambda, maxtheta, maxtime  I]  =  simp  le_reentry_reference(LDVALUE, altitude); 
global  edO  el  beta  betam  sigma  lambda  maxtime  1 
global  target 

global  hO  thetaO  phiO  vO  gamO  psiO  maxtheta  maxphi 
global  GO  SCALE  RHOO  MASS  LR  AREA  RE 


%  Times 
tstep  =  1 ;  %sec 
tfinal  =  4000;  %sec 


%  Earth  and  atmospherie  parameters 
RE  =  6378000;  GO  =  9.8; 

RHOO  =  1.752;  SCALE  =  6700; 


%  Vehiele  parameters 
cd0=  l;el  =  LDVALUE; 

beta  =  (MASS*G0)/(ed0*FR_AREA);  betam  =  beta/GO; 
sigma  =  sigmaopt(el,ed0); 

%  Initial  Conditions 

hO  =  altitude*  1000;  thetaO  =  0;  phiO  =  0; 
vO  =  sqrt(RE*G0)  -  130;  gamO  =  0;  psiO  =  0; 

%  Generate  referenee  trajeetory  with  sigma  optimal 

sigma  =  sigmaopt(ol,od0); 

yO  =  [hO  thetaO  phiO  vO  gamO  psiO]; 

options  =  odeset('RelTor,le-08,AbsTor,le-10*ones(l,6)); 

[t,y]  =  ODE45('simple_reentry_eoms',0:tstep:tfinal,yO, options); 

for  i  =  1  Tfinal; 
ify(U)>0 
p(i,l)  =  y(i,l); 
p(i,2)  =  y(i,2); 
p(i,3)  =  y(i,3); 
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p(i,4)  =  y(i,4); 
p(i,5)  =  y(i,5); 
p(i,6)  =  y(i,6); 
else  break 
end 
end 

[maxphi,maxtimel]  =  max(p(:,3)); 
lambda  =  maxphi; 
maxphi  =  maxphi*(  180/pi); 
maxtheta  =  p(maxtimel,2)*(180/pi); 
target  =  maxtheta  +  maxphi; 
missvalue  =  target  -  maxtheta; 

%  Now  try  to  aehieve  theta  =  target  by  using  differing  sigma  values  and 
%  roll  reversals.  If  the  SMV  is  short,  deerease  sigma;  if  long,  inerease 
%  sigma. 

sigma  =  sigmaopt(ol,odO); 

q  =  p; 

while  missvalue  >  .  1 ; 
clear  q; 

yO  =  [hO  thetaO  phiO  vO  gamO  psiO]; 

options  =  odeset('RelTol',le-08,'AbsTol',le-10*ones(l,6)); 

[t,y]  =  ODE45('reentry2eoms_nohft',0:tstep:tfinal,yO, options); 
for  i  =  1  :tfinal; 
ify(U)>0 
q(i,l)  =  y(i,l); 
q(i,2)  =  y(i,2); 
q(i,3)  =  y(i,3); 
q(i,4)  =  y(i,4); 
q(i,5)  =  y(i,5); 
q(i,6)  =  y(i,6); 
else  break 
end 
end 

[maxtheta2,maxtime2]  =  max(q(:,2)); 
maxtheta2  =  maxtheta2*(  180/pi); 
missvalue  =  target  -  maxtheta2; 
sigmacorr  =  missvalue; 
sigma  =  sigma  -  sigmacorr*(pi/180); 
end 

maxtime  1 
maxtime2 

timediff  =  abs(maxtimel  -  maxtime2) 
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%  Pass  maximum  time  back  to  main  program  for  use  in  coverage  statistics 
if  maxtime  1  >  maxtime2; 

downrange_time  =  (maxtheta*(pi/180))/m_motion; 
else 

downrangetime  =  maxtime2; 
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  REENTRY  EQUATIONS  OF  MOTION 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function  [ydot]  =  simple_reentry_eoms(t,y) 


global  cdO  cl  betam  sigma 

global  RE  GO  SCALE  RHOO 

global  hO  thetaO  phiO  vO  gamO  psiO  turnO 


if  y(6,l)  >=  psiO  +  turnO; 

sigma  =  0; 
end 


%dh/dt 

ydot(l,l)  =  -y(4)*sin(y(5)); 

%dtheta/dt 

ydot(2,l)  =  (y(4)*cos(y(5))*cos(y(6)))/(cos(y(3))*(RE  +  y(l))); 

%dphi/dt 

ydot(3,l)  =  (y(4)*cos(y(5))*sin(y(6)))/(RE  +  y(l)); 

%dv/dt 

ydot(4,l)  =  -(RHOO*exp(-y(l)/SCALE)/(2*betam))*y(4)^2  +  G0*sin(y(5)); 
%dgamma/dt 

ydot(5,l)  =  -(l/y(4))*((y(4)^2/(RE  +  y(l)))*cos(y(5))  +  ((RHOO* 
exp(-y(l)/SCALE)/(2*betam))*y(4)^2*(cEcdO))*cos(sigma)  -  G0*cos(y(5))); 

%dpsi/dt 

ydot(6,l)  =  (l/(y(4)*cos(y(5))))*((y(4)^2/(RE  +  y(l)))*(cos(y(5))^2)*sin(y(6))*tan(y(3)) 
+  ((RHOO*exp(-y(l)/SCALE)/(2*betam))*y(4)^2*(cl/cdO))*sin(sigma)); 
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Earth  Grid 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  THIS  PROGRAM  CREATES  A  GRID  OF  POINTS  SPACED  EQUALEY  AROUND 
%  THE  EARTH.  IT  ACCOUNTS  FOR  BUNCHING  AT  THE  POLES  AS  WELL  AS 
%  RANDOMIZING  THE  START  POSITION  OF  THE  FIRST  POINT  ALONG  EACH 
%  LINE  OF  LATITUDE. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


global  tol  num  lat  long  earth  grid  grid  tolerance  grid_points 
global  xl  yl  zl  max  lat 

%  User  sets  the  spacing  between  grid  points,  entered  in  degrees, 
tol  =  grid_tolerance*(pi/180); 

%  Initialize  to  zero  degrees  in  both  lat  and  long. 

%  Also  introduce  randomness  into  longitude  start  so  prime  meridian  doesn't 
%  get  too  much  weight, 
lat  =  0; 

long  =  0  +  rand*tol; 
num  =  0; 
a  =  1; 

%  Outer  loop  increments  latitude  and  keeps  longitude  starting  point  random, 
while  lat  <=  max  lat; 

%  Inner  loop  increments  longitude 
while  long  <=  2*pi; 

earth_grid(a,l)  =  lat*(180/pi); 
earth_grid(a,2)  =  long*(  180/pi); 
earth_grid(a,3)  =  0; 

%  This  block  takes  care  of  Southern  Hemisphere  by  mirroring  points 
%  when  latitude  is  not  zero, 
if  lat  >  0 
a  =  a  +  1; 

earth_grid(a,l)  =  -lat*(  180/pi); 
earth_grid(a,2)  =  long*(  180/pi); 
earth_grid(a,3)  =  0; 
end 

num  =  (2*pi/tol)*cos(lat); 
long  =  long  +  (2*pi/num); 
if  long  <=  2*pi 
a  =  a  +  1; 
end 
end 

long  =  0  +  rand*tol; 
lat  =  lat  +  tol; 
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a  =  a  +  1; 
end 

[grid_points,r]  =  size(earth_grid); 
gridindx  =  [l;grid_points]; 

current_lat  =  pi/2  -  (earth_grid(grid_indx,l).*(pi/180)); 
current_long  =  earth_grid(grid_indx,2).*(pi/180); 
xl  =  cos(current_long).*sin(current_lat); 
yl  =  sin(current_long).*sin(current_lat); 
zl  =  cos(current_lat); 


Constellation  Function 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  CREATES  A  NOMINAL  CONSTELLATION  BASED  ON  INPUTS  FROM  THE 
%  REENTRY  SIMULATION. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  [constellation, num_planes,num_sats]  = 

constellation_func2(lambda,maxtimel,inc,m_motion,n) 


global  constellation  earth  coverage  TTT  w 


%  Determine  number  of  SMVs  for  this  value  of  lambda 
smv_per_plane  =  ceil(2*pi/mod(m_motion*(TTT  -  maxtimel),2*pi)); 
if  inc  +  lambda  >=  pi/2; 

num_planes  =  ceil(pi/w); 
else 

num_planes  =  ceil((pi  +  n)/w); 
end 

num_sats  =  num_planes*smv_per_plane; 
raaninit  =  0; 
arglat  =  0; 
count  =  1 ; 

u_incr  =  2*pi/smv_per_plane; 
raanincr  =  w; 

constellation  =  zeros(num_sats,3); 
while  count  <  num  sats; 

for  count2  =  countxount  +  smv_per_plane  -  1; 
constellation(count2,l)  =  inc; 
constellation(count2,2)  =  raaninit; 
constellation(count2,3)  =  arglat; 
arg_lat  =  mod(arg_lat  +  u_incr,2*pi); 
end 
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count  =  count2  +  1 ; 
raaninit  =  raaninit  +  raan_incr; 
arglat  =  0; 
end 


Genetic  Algorithm 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  THIS  IS  THE  MAIN  PROGRAM  IN  THE  GENETIC  SEARCH.  IT  HAS  EOUR 
%  PARTS: 

%  I)  OBTAIN  NOMINAE  REENTRY  PARAMETERS  BASED  ON  SMV  INPUTS 
%  2)  BUIED  AN  INITIAL  POPULATION  OE  CONSTELLATIONS 
%  3)  GA  SEARCH  TO  FIND  OPTIMAL  CONSTELLATION 
%  4)  DISPLAY  RESULTS  OF  TWO  BEST  AND  TWO  WORST  CONSTELLATIONS 
%  GRAPHICALLY 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


clear  all;  paek; 

global  MASS  LDVALUE  FR_AREA  altitude  earth  grid  lambda  TTT  big  array  bigrows 
global  grid  toleranee  sma  ecc  ine  m  motion  MU  num  sats  fitness  A  fitfun  h  k 
global  hO  thetaO  phiO  vO  gamO  psiO  RE  GO  RHOO  SCALE  max  lat  maxgen  slope  yint 
global  planesmod  w 


%  Inputs  for  the  SMV  and  the  orbit,  as  well  as  the  tolerances  for  the  Earth 


%  grid  and  latitude  limits 
MASS  =  1000;  %  kg 

LDVALUE  =  .7;  %  lift/drag 

FR_AREA=10;  %  m^2 

deltav  =  0;  %  m/s 

altitude  =  500;  %  km 

max_lat  =  65*(pi/180);  %  deg 

grid_tolerance  =  5;  %  deg 

MU  =  3.986e5;  %  constant 

TTT  =  5400;  %  sec 

%  Nominal  orbital  elements 

sma  =  6378  +  altitude;  %  km 

ecc  =  0;  %  no  units 


%  Calculate  orbital  period 
m_motion  =  sqrt(MU/sma^3);  %  rad/s 
period  =  2*pi/m_motion;  %  see 


%  Get  displacement  distanee  tfom  referenee  reentry  profile 
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[lambda, maxtheta, maxtime  1  ]  =  simp  le_reentry_reference(LDVALUE, altitude); 

%  Set  up  and  calculate  Earth  grid 
earthgridpoints 

%  Generate  the  full  constellation  array  for  all  possible  SMVs 
binsize  =  20; 

big_array  =  zeros(2^binsize,4); 
incint  =  .05; 
raanint  =  .05; 
arglat_mt  =  .05; 

values  =  const_perms(inc_int,raan_int,arglat_int); 
values(:,4)  =  1; 

[valrows,valcols]  =  size(values); 
big_array(l:valrows,l:4)  =  values; 

[bigrows,bigcols]  =  size(big_array); 

%  Setup  max  generations  and  fitness  function  exponent 

maxgen  =  5000; 
maxexp  =  20; 
h=  l;k  =  3; 
fitfun  =  3; 

%  Einear 
if  fitfun  ==  1 ; 

slope  =  (maxexp  -  k)/(maxgen  -  h); 
yint  =  1  -  slope; 
end 

%  Right  parabola 
if  fitfun  ==  2; 

A  =  (maxgen  -  h)/((maxexp  -  k)^2); 
end 

%  Lip  parabola 
if  fitfun  ==  3; 

A  =  (maxexp  -  k)/((maxgen  -  h)^2); 
end 

%  Initialize  Genetic  Search 

seed  =  601387; 
gs_init(seed); 

rand('state',  91403) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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%  DEFINE  THE  INITIAL  POPULATION  BASED  ON  MAX  LAT.  INCLUDE 
%  HIGHER%  INCLINATIONS  EVENLY  SPACED  FROM  MAX  EAT  TO  90DEG. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  Define  the  range  of  inelinations  the  eonstellation  funetion  will  eonsider 
memcount  =  1; 

%  Define  the  baseline  SOC  eonstellation 
oount2  =  1 ; 

for  ine  =  [lambda;. 002:pi/2  -  lambda]; 
w  =  2*asin(sin(lambda)/sin(ino)); 
for  n  =  [0;.002:pi]; 

d  =  atan(sin(n/2)*tan(ino)); 

X  =  asin((sin(lambda)*oos(d))/oos(inc)); 
if  d  +  X  >=  max  lat; 
nplanes  =  (pi  +  n)/w; 
if  nplanes  >=  1; 
planes3(eount2)  =  (pi  +  n)/w; 
i3(eount2)  =  ine; 
n3(eount2)  =  n; 
eount2  =  eount2  +  1 ; 
end 
break 
end 
end 
end 

[planes_mod,p_indx]  =  min(planes3); 
m_ine  =  i3(p_indx); 
m_n  =  n3(p_indx); 

[fulll  ,mod_l  ,full_sats,mod_sats]  = 

eonstellation_func_mod(lambda, maxtime  1  ,m_inc,m_motion,m_n); 

%  Assign  index  numbers  tfom  big  array  to  eaeh  CAV  in  eaeh  eonstellation 
%  and  pad  the  remainder  of  the  ehromosome  with  zeros 
base_veot  =  zeros(l,100); 

[full_lrows,full_lools]  =  size(full_l); 
for  aa  =  Efull  lrows; 

aarow  =  find(big_array(:,l)  <=  full_l(aa,l)  +  ino_int/2  &  big_array(:,l)  >= 
full_l(aa,l)  -  ino_int/2  &  big_array(:,2)  <=  full_l(aa,2)  +  raan_int/2  & 
big_array(:,2)  >=  full_l(aa,2)  -  raan_int/2  &  big_array(:,3)  <=  full_l(aa,3)  + 
arglat_int/2  &  big_array(:,3)  >=  full_l(aa,3)  -  arglat_int/2); 
base_veot(aa)  =  aarow; 
end 

fulll  ehrom  =  gs_blank(reshape(dee2bin(base_vect,binsize)',  1 , 1 00*binsize)); 
mem  id  =  ['mem  id'  num2str(mem_eount)]; 
memid  =  gs_new('Popr,full_l  ehrom) 
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memcount  =  memcount  +  1; 

base_vect  =  zeros(l,100); 

[mod_lrows,mod_lcols]  =  size(mod_l); 
for  aa  =  l:mod_lrows; 

aarow  =  find(big_array(:,l)  <=  mod_l(aa,l)  +  inc_int/2  &  big_array(:,l)  >= 
mod_l(aa,l)  -  inc_int/2  &  big_array(:,2)  <=  mod_l(aa,2)  +  raan_int/2  & 
big_array(:,2)  >=  mod_l(aa,2)  -  raan_int/2  &  big_array(:,3)  <=  mod_l(aa,3)  + 
arglat_int/2  &  big_array(:,3)  >=  mod_l(aa,3)  -  arglat_int/2); 
base_vect(aa)  =  aarow; 
end 

mod_lchrom  =  gs_blank(reshape(dec2bin(base_vect,binsize)',l,100*binsize)); 
mem  id  =  ['mem  id'  num2str(mem_count)]; 
memid  =  gs_new('Popr,mod_lchrom) 
memcount  =  memcount  +  1; 

%  Now  add  constellations  from  a  range  of  inclinations 
inc_fidelity  =  (i3(end)  -  i3(l))/length(i3); 
incincrement  =  20; 

inc_spacing  =  (pi/2  -  min(i3))/inc_increment; 
inccounter  =  0; 
nextindx  =  []; 

while  inc  counter  <  inc  increment; 
while  isempty(next_indx)  ==  1 ; 

nextindx  =  find(i3(l  ,0  <=  (i3(l)  +  inc_spacing*inc_counter)  +  inc_fidehty  &... 

i3(l,:)  >=  (i3(l)  +  inc_spacing*inc_counter)  -  inc  fidelity); 
incfidelity  =  incfidelity  +  .001; 
end 

nextinc  =  i3(next_indx(l)); 
planes_mod  =  planes3(next_indx(l)); 
nextn  =  n3(next_indx(l)); 

[fulll  ,mod_l  ,full_sats,mod_sats]  = 
constellation_func_mod(lambda, maxtime  l,next_inc,m_motion,next_n); 

%  Assign  index  numbers  from  big  array  to  each  CAV  in  each  constellation 
%  and  pad  the  remainder  of  the  chromosome  with  zeros 
basevect  =  zeros(l,100); 

[full_lrows,full_lcols]  =  size(full_l); 
for  aa  =  l:full_lrows; 

aarow  =  find(big_array(:,l)  <=  full_l(aa,l)  +  inc_int/2  &  big_array(:,l)  >= 
full_l(aa,l)  -  inc_int/2  &  big_array(:,2)  <=  full_l(aa,2)  +  raan_int/2  & 
big_array(:,2)  >=  full_l(aa,2)  -  raan_int/2  &  big_array(:,3)  <=  full_l(aa,3)  + 
arglat_int/2  &  big_array(:,3)  >=  full_l(aa,3)  -  arglat_int/2); 
base_vect(aa)  =  aarow; 
end 
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full_l  chrom  =  gs_blaiik(reshape(dec2bin(base_vect,binsize)',  1 , 1 00*binsize)); 
mem_id  =  ['mem  id'  num2str(mem_count)]; 
memid  =  gs_new('Popr,full_l  chrom) 
memcount  =  memcount  +  1 ; 

basevect  =  zeros(l,100); 

[mod  i  rows, mod  i  cols]  =  sizc(mod_l); 
for  aa  =  l:mod_lrows; 

aarow  =  find(big_array(:,l)  <=  mod_l(aa,l)  +  inc_int/2  &  big_array(:,l)  >= 
mod_l(aa,l)  -  inc_int/2  &  big_array(:,2)  <=  mod_l(aa,2)  +  raan_int/2  & 
big_array(:,2)  >=  mod_l(aa,2)  -  raan_int/2  &  big_array(:,3)  <=  mod_l(aa,3)  + 
arglat_int/2  &  big_array(:,3)  >=  mod_l(aa,3)  -  arglat_int/2); 
basc_vcct(aa)  =  aarow; 
end 

mod_  1  chrom  =  gs_blarLk(reshape(dcc2bin(base_vect,binsize)',  1 , 1 00*binsize)); 
mem_id  =  ['mem  id'  num2str(mem_count)]; 
memid  =  gs_new('Popr,mod_l  chrom) 
memcount  =  memcount  +  1 ; 

%  Increment  latitude  loop  and  return  to  top  of  section 
inccounter  =  inc_counter  +  1; 
inc_fidelity  =  (i3(end)  -  i3(l))/length(i3); 
nextindx  =  []; 
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  CREATE  RANDOM  MEMBERS  FROM  BIG  ARRAY 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


seed_array  =  randperm(bigrows); 
linecount  =  1 ; 

for  jj  =  mem_count  +  l:mem_count  +101; 
randlength  =  floor(rand*100); 
nextehrom  =  zeros(  1 , 1 00); 

next_chrom(l,l;l  +  randlength)  =  seed_array(l, linecount:  linecount  +  randlength); 
nextehrom  =  gs_blank(reshape(dec2bin(next_chrom,binsize)',  1 , 1 00*binsize)); 
linecount  =  linecount  +  randlength; 

IDstr  =  ['mem_id'  num2str(jj)]; 

IDstr  =  gs_new('Popr,next_chrom); 
end 

%  Throw  in  some  "all  zeros"  chromosomes  for  variety... 
zerochrom  =  gs_blarLk(num2str(zeros(  1,2000))); 
for  kk  =  1:100; 

IDstr  =  ['mem  id'  num2str(jj  +  kk)]; 
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IDstr  =  gs_new('Popr,zero_chrom); 
end 

%  Add  some  "all  ones"  ehromosomes  for  even  more  variety... 
oneschrom  =  gs_blarLk(num2str(ones(  1,2000))); 
for  11=  1:100; 

IDstr  =  ['mem_id'  num2str(jj  +  kk  +  11)]; 

IDstr  =  gs_new('Popr,ones_ohrom); 
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  MAIN  GENETIC  MANIPULATION  LOOP 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  Launch  interrupt  buttons 

gs_open_cbox; 

%  Evaluate  initial  population 

memcount  =  gs_popsize('Popr); 

for  id  =  Lmem  count 

chrl  =  gs_get('Popr,  id); 

fitness  =  main_lltness_func_bin(chrl,3) 

mem  id  =  gs_set_fit('Popr,  id,  fitness); 


%  Check  for  suspend  and  break  signals 

gs_break; 

gs_suspend; 


end 

%  Lind  best  member  in  the  initial  population 

disp('The  best  member  in  the  initial  population  is'); 
memids  =  gs_sel_lofit('Popr); 
chrl  =  gs_get('Popr,  mem_ids(l)) 
disp('Its  fitness  is') 

fitness  =  gs_get_fit('Popr,  mem_ids(l)) 
best_lit_start  =  fitness; 

%Save  the  best  initial  constellation  for  comparison  later... 
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[ee]  =  gs_sel_lofit('Popr); 

[fitness, earth_coverage]  =  main_fitness_func_bin(gs_get('Popr,ee(l)),maxgen) 
aa  =  bin2dec(reshape(gs_unblaiik(gs_get('Pop  1  ',ee(2))),20, 1 00)')'; 
aa  ~  aa'; 

lopopinitial  =  big_array(aa(find(aa), :),:); 
displayinitial  =  lopop_initial(:,l;3)*(180/pi) 
pause; 

%  Cheek  for  break  signal 
gs_break; 

%  Genetie  Seareh  Loop 

for  gen  =  Lmaxgen 
gen 

%  Trim  population  if  over  500  members 

memeount  =  gs_popsize('Popr); 
if  mem  eount  >500 

mem_ids  =  gs_selr_hifit('Popr); 
gs_del('Pop  1  ',mem_ids(  1 )); 
if  mem_oount-l  >  500 

gs_del('Pop  1  ',mem_ids(2)); 
end 
end 

%  Select  genetic  operation 

op  name  =  gs_sel_op({'mutbin',  'xovr2'},  [0.200000,  0.800000]); 

%  Select  members  for  genetic  operation 

switch  char(op_name) 
case  'mutbin' 

memids  =  gs_selr_lofit('Popr); 
case  'xovr2' 

mem_ids  =  gssehCPopl'); 
end 

%  Implement  genetic  operation 

off  ids  =  gs  opCPopT,  op  name,  mem_ids(l),  mem_ids(2),  0.200000); 
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%  Evaluate  the  fitness  of  the  offspring 

for  off  =  l:length(off_ids) 

ehrl  =  gs_get('Popr,off_ids(off)); 
fitness  =  main_fitness_func_bin(ehrl,gen) 
mem  id  =  gs_set_fit('Popr,  off_ids(off),  fitness); 


end 

%  Check  for  suspend  and  break  signals 

gs_break; 

gs_suspend; 


end 

%  Find  the  best  individual 

disp('The  member  with  the  best  fitness  is') 
memids  =  gs_sel_lofit('Popl'); 
ehrl  =  gs_get('Popl',  mem_ids(l)) 
disp('Its  fitness  is'); 

fitness  =  gs_get_fit('Popl',  mem_ids(l)) 
best_fit_end  =  fitness; 

%  Close  interrupt  buttons 

gsclosecbox 

%  Show  how  much  GA  was  able  to  improve  over  the  intial  population 
improvement  =  best  fit  start  -  best  fit  end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  THIS  SECTION  PRINTS  A  SINGEE  CONSTEEEATION  ALONG  WITH  THE 
%  ORBITAL  ELEMENTS  OF  THE  CAVS  WITHIN  IT. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


[ee]  =  gs_sel_lofit('Popr); 

[fitness, earth_coverage]  =  main_fitness_func_bin(gs_get('Popr,ee(I)),maxgen) 
aa  =  bin2dec(reshape(gs_unblank(gs_get('Pop  I  ',ee(  I  ))),20, 1 00)')'; 
aa  ~  aa'; 

lopop  =  big_array(aa(find(aa), :),:); 

[losats, nothing]  =  size(lopop); 
const_elements_final  =  lopop(:,I:3).*(I80/pi) 
for  cur  sat  =  ITosats; 
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%  Get  orbital  elements  for  next  SMV 
incnow  =  lopop(cur_sat,l); 
raannow  =  lopop(eur_sat,2); 
arglatnow  =  lopop(cur_sat,3); 

%  Calculate  values  outside  loop  to  improve  speed 
cosine  =  cos(inc_now); 
sininc  =  sin(inc_now); 
count  =  1 ; 

%  Increment  footprint  until  90  minute  limit 

for  tt  =  arglat  now  +  maxtheta*(pi/180):.05:arglat_now  +  m_motion*(TTT  - 
maxtimel)+  maxtheta*(pi/180); 

%  Update  the  SSP  for  this  time  step  based  on  SMV's  orbital  elements 
SSPlat  =  asin(sin_inc*sin(tt)); 

SSPlong  =  mod(atan(cos_inc*tan(tt))  +  raan_now,2*pi); 
if  mod(tt,2*pi)  >  pi/2  &&  mod(tt,2*pi)  <  3*pi/2; 

SSP  long  =  mod(SSP_long  -  pi,2*pi); 
end 

SSP(count,l)  =  SSP_long*(  180/pi); 

SSP(count,2)  =  SSP_lat*(  180/pi); 
count  =  count  +  1 ; 
end 

hold  on; 

figure(l);  plot(SSP(:,l),SSP(:,2),'k.') 

strl  =  ['SMVs=  '  num2str(losats)  Coverage  =  '  num2str(earth_coverage)  ... 

L/D=  '  nuni2str(LDVALUE)  \phi=  '  num2str(max_lat*  180/pi)  ... 
inc=  '  num2str(inc_now*(  180/pi))  '\circ']; 
axis([0  360  -90  90]);  grid  off;  box  on; 
xlabel('Longitude');  ylabel('Latitude'); 

text(180,-85,strl,'HorizontalAhgnmentVcenterVBackgroundColorVwVEdgeColorVk'); 

end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  THIS  SECTION  PRINTS  THE  BEST  INITIAL  CONSTELLATION  ALONG  WITH 
%  THE  ORBITAL  ELEMENTS  OE  THE  CAVS  WITHIN  IT. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


[losats, nothing]  =  size(lopop_initial); 
constelementsinitial  =  lopop_initial(:,I:3).*(I80/pi) 
for  cur  sat  =  ITosats; 

%  Get  orbital  elements  for  next  SMV 
incnow  =  lopop_initial(cur_sat,l); 
raannow  =  lopop_initial(cur_sat,2); 
arglatnow  =  lopop_initial(cur_sat,3); 

%  Calculate  values  outside  loop  to  improve  speed 
cosine  =  cos(inc_now); 
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sininc  =  sin(inc_now); 
count  =  1 ; 

%  Increment  footprint  until  90  minute  limit 

for  tt  =  arglat_now  +  maxtheta*(pi/180):.05;arglat_now  +  m_motion*(TTT  - 
maxtimel)+  maxtheta*(pi/180); 

%  Update  the  SSP  for  this  time  step  based  on  SMV's  orbital  elements 
SSPlat  =  asin(sin_inc*sin(tt)); 

SSPlong  =  mod(atan(cos_ine*tan(tt))  +  raan_now,2*pi); 
if  mod(tt,2*pi)  >  pi/2  &&  mod(tt,2*pi)  <  3*pi/2; 

SSP  long  =  mod(SSP_long  -  pi,2*pi); 
end 

SSP(eount,l)  =  SSP_long*(  180/pi); 

SSP(count,2)  =  SSP_lat*(  180/pi); 
count  =  count  +  1 ; 
end 

hold  on; 

figure(2);  plot(SSP(:,l),SSP(;,2),'k.') 

str2  =  ['SMVs=  '  num2str(losats)  Coverage  =  '  num2str(earth_coverage)  ... 
L/D=  '  num2str(LDVALUE)  \phi=  '  num2str(max_lat*  180/pi)  ... 
inc=  '  num2str(ine_now*(  180/pi))  '\cire']; 
axis([0  360  -90  90]);  grid  off;  box  on; 
xlabel('Longitude');  ylabel('Latitude');  text(180,- 
85,str2,'HorizontalAlignmentVcenterVBackgroundColorVwVEdgeColorVk'); 


Constellation  Fitness  Function 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  THIS  FUNCTION  DETERMINES  THE  FITNESS  OF  A  CONSTEEEATION  BY 
%  COMBINING  THE  NUMBER  OF  CAVS  WITH  EARTH  COVERAGE. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  [fitness, earth  coverage]  =  main_fitness_func_bin(constellation,gen) 
global  earth  grid  lambda  grid_points  A  fitfun 
global  m  motion  TTT  big  array  slope  yint  h  k 
global  maxtheta  maxtime  1  xl  yl  zl  maxgen  bigrows 


%  Apply  appropriate  model  for  exponent 
%  Einear 
if  fitfun  ==  1 ; 

n  =  slope*gen  +  yint; 
end 

%  Right  parabola 
if  fitfun  ==  2; 
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n  =  sqrt((gen  -  h)/A)  +  k; 
end 

%  Up  parabola 
if  fitfun  ==  3; 

n  =  A*  (gen  -  h)^2  +  k; 
end 

%  Reshape  the  chromosome  string  into  a  matrix  representing  the 
%  constellation  to  be  evaluated 

vectl  =  birL2dec(reshape(gs_unblarLk(constellation),20,100)')'; 

if  sum(vectl,2)  ==  0; 
numsats  =  1 ; 
earthcoverage  =  0.1; 
fitness  =  num_sats/(earth_coverage^n); 
return 
end 

vectl  =  vectl'; 

chrl_const  =  big_array(vectl(find(vectl 
[numsats, nothing]  =  size(chrl_const); 
earth_grid(:,3)  =  0; 
for  chrom  count  =  Imum  sats; 

if  chrl_const(chrom_count,4)  ==  1; 

%  Get  orbital  elements  for  next  SMV 
incnow  =  chrl_const(chrom_count,l); 
raan_now  =  chrl_const(chrom_count,2); 
arglatnow  =  chrl_const(chrom_count,3); 

%  Calculate  values  outside  loop  to  improve  speed 
cos_inc  =  cos(inc_now); 
sin_inc  =  sin(inc_now); 
swath_test  =  cos(lambda); 

%  Set  stepsize  for  propagation  loop 
ttincr  =  lambda; 

%  Increment  footprint  until  90  minute  limit 

for  tt  =  arglat  now  +  maxtheta*(pi/180):tt_incr:arglat_now  +  m_motion*(TTT  - 
maxtimel)+  maxtheta*(pi/180); 

%  Update  the  SSP  for  this  time  step  based  on  SMV's  orbital  elements 
SSP_lat  =  pi/2  -  asin(sin_inc*sin(tt)); 

SSPlong  =  mod(atan(cos_inc*tan(tt))  +  raan_now,2*pi); 
if  mod(tt,2*pi)  >  pi/2  &&  mod(tt,2*pi)  <  3*pi/2; 

SSP  long  =  mod(SSP_long  -  pi,2*pi); 
end 

x2  =  cos(SSP_long)*sin(SSP_lat); 
y2  =  sin(SSP_long)*sin(SSP_lat); 


67 


z2  =  cos(SSP_lat); 

%  Perform  angular  distance  check  of  grid  point  and  update  grid  coverage 
testl  =  (xl*x2  +  yl*y2  +  zl*z2); 
indxl  =  find(testl  >=  swath  test); 

earth_grid(indxl,3)  =  earth_grid(indxl,3)  +  1; 
end 
end 
end 

%  Calculate  Earth  coverage 
coveragecounter  =  0; 
xtracov  =  0; 
for  cc  =  l;grid_points; 
if  earth_grid(cc,3)  >=  1; 

coveragecounter  =  coveragecounter  +  1; 
end 

xtra_cov  =  xtra_cov  +  earth_grid(cc,3); 
end 

earth_coverage  =  coverage_counter/grid_points; 
extra_coverage  =  xtra_cov/grid_points; 

%  Return  fitness  for  this  constellation 
if  earth_coverage  >=  .5; 

fitness  =  num_sats/(earth_coverage^n); 
else 

fitness  =  num_sats/(.l^n); 
end 
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